
# 加载所需库------
if (!require("pacman")) install.packages("pacman")
pkgs <- c("dplyr", "ggplot2", "readxl", "tidyr", "lubridate", "janitor", "stringr", "writexl",
          "psych", "magrittr", "sf", "spdep", "tmap", "ggspatial","RColorBrewer")
pacman::p_load(pkgs, character.only = TRUE)

# 1.读取数据------
# 1.1读取发病数数据------
data <- read_excel("./data/processed/cleaned_data9.xlsx") %>%
  filter(year >= 2004)
data$city_code %>% unique()
data$county_code %>% unique()

#读取中英文名字
source("./code/2.3 人口数估发病率.R")

# 1.2读取发病率数据（系统导出）------
rate <- readxl::read_excel("data/origin/发病率2004-2023年v1.xlsx", 
           sheet = "Report", col_types = c("text", 
                                           "text", "numeric", "numeric", "numeric", 
                                           "numeric", "numeric", "numeric", 
                                           "numeric", "numeric", "numeric", 
                                           "numeric", "numeric", "numeric", 
                                           "numeric", "numeric", "numeric", "numeric", 
                                           "numeric", "numeric", "numeric", "text")) %>% 
  filter(region != "county兵团")
warnings()
rate$地区编码 %>% str_sub(start = 1, end = 4) %>% unique()
rate$地区 %>% unique()
rate$region %>% unique();1+96+14

# 重命名和整理列名
colnames(rate)[3:22] <- gsub("年", "", colnames(rate)[3:22]) # 去掉年份列名中的“年”字

rate2 <- rate %>%
  ungroup() %>%
  pivot_longer(cols = -c(地区, 地区编码, region), names_to = "year", values_to = "incidence_rate") %>%
  mutate(
    incidence_rate = ifelse(incidence_rate == "-", "0", incidence_rate), # 替换"-"为"0"
    incidence_rate = as.numeric(incidence_rate), # 转换为数值型
    year = as.integer(year), # 确保年份为整数
    county_code = substr(地区编码, 1, 6), # 提取县区代码
    city_code = substr(地区编码, 1, 4)) # 提取城市代码

rate2$city_code %>% unique()

# 重新编码：6540变为6541；654203/654223
rate2 <- rate2 %>%
  mutate(
    city_code = case_when(
      city_code == "6541" ~ "6540",
      TRUE ~ city_code # 默认情况下保留原来的值
    ),
    county_code = case_when(
      county_code == "654203" ~ "654223",
      TRUE ~ county_code # 默认情况下保留原来的值
    ))
rate2$county_code <- sub("^6541", "6540", rate2$county_code)
rate2<- rate2 %>%  filter(region != "county兵团")
dim(rate2)

# 1.3读取发病率数据(人口数和发病数估计）------

rate2_estimate <- rate_year %>% 
  mutate(county_code = 地区编码,city_code = substr(county_code,1,4)) %>%
  dplyr::rename(incidence_rate = rate) %>%
  dplyr::select(地区,地区编码, region,year,incidence_rate,county_code,city_code)
table(rate2$region);table(rate2_estimate$region);14*19;96*19;1*19
dim(rate2);dim(rate2_estimate)


wb1 <- createWorkbook()  # 新建空工作簿（无任何 Sheet）
if (!"发病率_年" %in% names(wb1)) {addWorksheet(wb1, "发病率_年")}
writeData(wb1, sheet = "发病率_年", x = rate_year)


### 2.1 整个新疆地区2005-2023年年度发病数和发病率------
# 在两个图的scale_fill_gradient()中添加以下参数
scale_fill_gradient(
  low = "white", high = "#FB8072",
  limits = c(0, 150),  # 强制图例范围为0-150
  breaks = c(0, 50, 100, 150),  # 指定断点
  labels = c("0", "50", "100", "150"))  # 自定义标签（可选）

# rate2 <- rate2_estimate
psych::describe(rate2_estimate)
psych::describe(rate2)

rate2 <- rate2_estimate

rate2$city_code %>% unique()
rate2$county_code %>% unique() %in% "654203"

rate2$county_code <- recode(rate2$county_code, `652201` = '650502', `654203` = '654223')

#Mann-Kendall趋势检验
rate_xj <- rate2 %>% filter(region == "province")

# 安装并加载所需的包
# if (!require(trend)) {
#   install.packages("trend")}
# library(trend)

# 提取年份和发病率数据
dataxj <- rate_xj %>%
  select(year, incidence_rate) %>%
  arrange(year)  # 确保按年份排序

# 执行 Mann-Kendall 趋势检验,使用trend包
# mk_test <- trend::mk.test(dataxj$incidence_rate)
# print(mk_test)

# 计算年度发病数和发病率
yearly_cases <- data %>%
  filter(year >= 2005, year <= 2023) %>%
  count(year, name = "cases")

# yearly_incidence <- rate2 %>%用的旧数据2004年系统导出的数据
#   filter(地区 == "新疆") %>%
#   select(year, incidence_rate) %>%
#   left_join(yearly_cases, by = "year")

yearly_incidence <- rate_year_province %>% dplyr::select(year,incidence_rate = rate,cases = count)

writexl::write_xlsx(yearly_incidence, "./output/表格2.1_发病数和发病率_年度.xlsx")


### 绘图：整个新疆地区2005-2023年年度发病数和发病率（双Y轴图）

# 绘制双Y轴图
ggplot(yearly_incidence) +
  geom_col(aes(x = year, y = cases), fill = "#8DD3C7", alpha = 0.6) +
  geom_line(aes(x = year, y = incidence_rate * 1000 / 4), color = "#FB8072", size = 1) +
  geom_point(aes(x = year, y = incidence_rate * 1000 / 4), color = "#FB8072") +
  scale_y_continuous(
    name = "Case Count",
    sec.axis = sec_axis(~ . / 1000 * 4, name = "Incidence Rate (per 100,000)")
  ) +
  labs(title = "Annual Case Count and Incidence Rate in Xinjiang (2005-2023)", x = "Year") +
  theme_minimal() +
  theme(
    axis.title.y = element_text(color = "#8DD3C7"),
    axis.title.y.right = element_text(color = "#FB8072"),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )
# ggsave("./主图/图2.2 整个新疆地区2005-2023年年度发病数和发病率.png")

### 2.2 分区县，1月-12月各月份的发病数和发病率------

# 读取和修复地图
rate2_estimate$region %>% unique()
name_county <- rate2_estimate %>% filter(region == "county")%>%
  dplyr::select(county_code,地区) %>% unique() %>% left_join(nameCE_county)
name_city <- rate2_estimate %>% filter(region == "city")%>%
  dplyr::select(city_code,地区) %>% unique() %>% left_join(nameCE_city)

# 计算1月-12月各月份的发病数和发病率
monthly_cases_by_county <- data %>%
  filter(year >= 2004, year <= 2023) %>%
  group_by(county_code, month) %>%
  summarise(total_cases = n(), .groups = 'drop') %>%
  left_join(name_county, by = "county_code")
monthly_cases_by_county$county_code %>% unique()

monthly_cases_ratio <- monthly_cases_by_county %>%
  group_by(month) %>%
  summarise(total_monthly_cases = sum(total_cases), .groups = 'drop') %>%
  mutate(total_all_cases = sum(total_monthly_cases), ratio = total_monthly_cases / total_all_cases)

print(monthly_cases_ratio)

monthly_cases_ratio


# 计算各月份累计发病数和发病率
monthly_cases_ratio <- rate_month %>%
  filter(region == "province") %>%
  group_by(month) %>%
  summarise(
    total_monthly_cases = sum(count),  # 各月份累计发病数
    incidence_rate = (sum(count) / mean(value)) * 100000  # 各月份发病率
  ) %>%
  mutate(month = factor(month, levels = 1:12))  # 将月份转换为因子，确保顺序正确

# 1. 计算各月份累计发病数和发病率
monthly_cases_ratio <- rate_month %>%
  filter(region == "province") %>%
  group_by(month) %>%
  summarise(
    total_monthly_cases = sum(count),  # 各月份累计发病数
    incidence_rate = (sum(count) / mean(value)) * 100000  # 各月份发病率
  ) %>%
  mutate(month = factor(month, levels = 1:12))  # 将月份转换为因子，确保顺序正确

writexl::write_xlsx(monthly_cases_ratio, "./output/表格2.1_发病数和发病率_月份累积.xlsx")

# 2. 绘制双Y轴图
# 计算两个轴的最大值，确保次Y轴能够正确显示
max_cases <- max(monthly_cases_ratio$total_monthly_cases)
max_rate <- max(monthly_cases_ratio$incidence_rate)
scale_factor <- max_cases / max_rate

ggplot(monthly_cases_ratio) +
  geom_col(aes(x = month, y = total_monthly_cases), fill = "#8DD3C7", alpha = 0.6) +
  geom_point(aes(x = month, y = incidence_rate * scale_factor), color = "#FB8072") +
  geom_line(aes(x = month, y = incidence_rate * scale_factor), color = "#FB8072", size = 1,group = 1) +
  scale_y_continuous(
    name = "Case Count",
    sec.axis = sec_axis(~ . / scale_factor, name = "Incidence Rate (per 100,000)")
  ) +
  labs(title = "Monthly Case Count and Incidence Rate in Xinjiang (2005-2023)", x = "Month") +
  theme_minimal() +
  theme(
    axis.title.y = element_text(color = "#8DD3C7"),
    axis.title.y.right = element_text(color = "#FB8072"),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )


# 保存图表
# ggsave("./主图/图2.3 整个新疆地区2005-2023年各月发病数和发病率.png")


if (!"发病率_月" %in% names(wb1)) {addWorksheet(wb1, "发病率_月")}
writeData(wb1, sheet = "发病率_月", x = monthly_cases_ratio)



### 2.3 分地市1月-12月各月份的发病率热力图------

# 确保 data 数据框中的 city_code 和 month 列存在
if (!all(c("city_code", "month") %in% colnames(data))) {
  stop("Columns 'city_code' or 'month' do not exist in the data frame. Please check the column names.")
}

# 2.4.1 计算每个月各城市的累计病例数
monthly_cases_by_city <- data %>%
  filter(year >= 2004, year <= 2023) %>%
  group_by(city_code, month) %>% #累计
  summarise(total_cases = n(), .groups = 'drop') %>%
  complete(city_code, month, fill = list(total_cases = 0)) %>%
  arrange(desc(total_cases))%>%
  left_join(name_city, by = "city_code")
monthly_cases_by_city$city_code %>% unique()

writexl::write_xlsx(monthly_cases_by_city, "./output/表格2.4.1_发病数和发病率_城市.xlsx")

if (!"发病率_城市_年月" %in% names(wb1)) {addWorksheet(wb1, "发病率_城市_年月")}
writeData(wb1, sheet = "发病率_城市_年月", x = monthly_cases_by_city)

# 提取 city_code 的排序顺序
city_code_order <- monthly_cases_by_city$city_code[!duplicated(monthly_cases_by_city$city_code)]

# 将地区变量转换为因子，并指定与 city_code 一致的顺序
# monthly_cases_by_city$地区 <- factor(
#   monthly_cases_by_city$地区,
#   levels = unique(monthly_cases_by_city$地区[match(city_code_order, monthly_cases_by_city$city_code)])
# )

# 绘制热力图（月发病数）
ggplot(monthly_cases_by_city, aes(x = factor(month), y = nameE, fill = total_cases)) +
  geom_tile() +
  scale_fill_gradient(low = "white", high = "#FB8072") +  # 使用Set3中的红色
  labs(
    # title = "Monthly Cumulative Cases Heatmap by City",
    x = "Month",
    y = "City",
    fill = "Total Cases"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )
# ggsave("./figs/Figure_2.4.1_Monthly_Incidence_Heatmap（月发病数）.png", width = 10, height = 6, dpi = 300)


# 绘制热力图（月发病率）
rate_month_province <- rate_month %>% filter(region =="province") %>%
  # mutate(city_code = substr(地区编码,1,4)) %>%
  group_by(地区,  month) %>%
  summarise(total_cases = sum(count), pop = mean(value)) %>%
  mutate(incidence  = total_cases/pop * 100000) %>%   ungroup()

rate_month_city <- rate_month %>% filter(region =="city") %>%
  mutate(city_code = substr(地区编码,1,4)) %>%
  group_by(city_code, 地区,  month) %>%
  summarise(total_cases = sum(count), pop = mean(value)) %>%
  mutate(incidence  = total_cases/pop * 100000) %>%   ungroup() %>%
  left_join(nameCE_city)

ggplot(rate_month_city, aes(x = factor(month), y = nameE, fill = incidence)) +
  geom_tile() +
  scale_fill_gradient(
    low = "white", high = "#FB8072",
    limits = c(0, 150),
    breaks = c(0, 50, 100, 150),
    labels = c("0", "50", "100", "150"),
    oob = scales::squish,
    name = "Incidence (/100,000)"# 直接通过梯度标度命名更规范
  ) +
  labs(
    # title = "Monthly Cumulative Incidence Heatmap by City",
    x = "Month",
    y = "Region",
    fill = "Incidence (/100 000)" 
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.title = element_text(size = 12, face = "bold")  # 增强标题可读性[1,3](@ref)
  )


# 保存图片
ggsave("./figs/Figure_2.4.1_Monthly_Incidence_Heatmap（月发病率）.png", width = 10, height = 6, dpi = 300)
ggsave("./主图/月发病率 热力图.png", width = 10, height = 6, dpi = 300)


# 2.4.1 计算各年各城市的累计病例数
yearly_cases_by_city <- data %>%
  filter(year >= 2004, year <= 2023) %>%
  group_by(city_code, year) %>%
  summarise(total_cases = n(), .groups = 'drop') %>%
  complete(city_code, year, fill = list(total_cases = 0)) %>%
  arrange(year)  %>% left_join(name_city, by = "city_code") %>%
  left_join(nameCE_city)
writexl::write_xlsx(yearly_cases_by_city, "./output/表格2.4.1_发病数和发病率_各城市各月份.xlsx")

if (!"发病率_城市_年" %in% names(wb1)) {addWorksheet(wb1, "发病率_城市_年")}
writeData(wb1, sheet = "发病率_城市_年", x = yearly_cases_by_city)

# 提取 city_code 的排序顺序
city_code_order <- yearly_cases_by_city$city_code[!duplicated(yearly_cases_by_city$city_code)]

# 将地区变量转换为因子，并指定与 city_code 一致的顺序
yearly_cases_by_city$地区 <- factor(
  yearly_cases_by_city$地区,
  levels = unique(yearly_cases_by_city$地区[match(city_code_order, yearly_cases_by_city$city_code)])) 

# 绘制热力图（年发病数）
ggplot(yearly_cases_by_city, aes(x = factor(year), y = nameE, fill = total_cases)) +
  geom_tile() +
  scale_fill_gradient(low = "white", high = "#FB8072") +  # 使用Set3中的红色
  labs(
    title = "Annual Cumulative Cases Heatmap by City",
    x = "Year",
    y = "Region",
    fill = "Total Cases") +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
    plot.title = element_text(hjust = 0.5, face = "bold"))

# 保存图片
# ggsave("./figs/Figure_Annual_Incidence_Heatmap(年发病数).png", width = 10, height = 6, dpi = 300)

# 绘制热力图（年发病率）
rate_year_province <- rate_year %>%
  ungroup() %>%filter(region == "province") %>% 
  mutate(rate = as.numeric(rate)) 

rate_year_city <- rate_year %>%
  ungroup() %>%filter(region == "city") %>% 
  mutate(rate = as.numeric(rate)) %>%
  left_join(name_city)

ggplot(rate_year_city, aes(x = factor(year), y = nameE, fill = rate)) +
  geom_tile() +
  scale_fill_gradient(
    low = "white", high = "#FB8072",
    limits = c(0, 150),
    breaks = c(0, 50, 100, 150),
    labels = c("0", "50", "100", "150"),
    oob = scales::squish,
    name = "Incidence (/100,000)")+ # 直接通过梯度标度命名更规范
  labs(
    title = "Annual Cumulative Incidence Heatmap by Region",
    x = "Year",
    y = "Region",
    fill ="Incidence (/100 000)") +
    theme_minimal() +
    theme(
      axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
      plot.title = element_text(hjust = 0.5, face = "bold"),
      legend.title = element_text(size = 12, face = "bold")  # 增强标题可读性[1,3](@ref)
    )

# 保存图片
ggsave("./figs/Figure_Annual_Incidence_Heatmap(年发病率).png", width = 10, height = 6, dpi = 300)
ggsave("./主图/年发病率 热力图.png", width = 10, height = 6, dpi = 300)


# 3、地区分布：分区县发病数、发病率、空间自相关_-----

# 3.1. 2005-2023年每年的发病数和发病率（作图用）------
# yearly_cases_county <- data %>%
#   filter(year >= 2004, year <= 2023) %>%
#   group_by(year, county_code) %>% 
#   summarise(cases = n(), .groups = 'drop') %>%
#   complete(year, county_code, fill = list(cases = 0)) %>% ungroup() %>% 
#   mutate(city_code = substr(county_code, 1, 4))
# names(yearly_cases_county)
yearly_cases_county <- rate_year_county %>% 
  rename(county_code = 地区编码) %>%
  mutate(city_code = substr(county_code, 1, 4)) %>%
  dplyr::select(year, county_code, cases = count, city_code) %>% 
  arrange(year, county_code)


#排除兵团
rate2_ex_bt <- rate2 %>% filter(!city_code %in% c("6590", "6500")) %>% 
  filter(region == "county") 

rate2_ex_bt$city_code %>% unique();yearly_cases_county$city_code %>% unique()
rate2_ex_bt$county_code %>% unique()
yearly_cases_county$county_code %>% unique()

yearly_incidence_county <- rate2_ex_bt %>% 
  filter(region == "county") %>%
  left_join(yearly_cases_county, by = c("year", "county_code","city_code")) %>%
  select(-地区编码) %>% arrange(year)
yearly_incidence_county$county_code %>% unique()
yearly_incidence_county$city_code %>% unique()

writexl::write_xlsx(yearly_incidence_county, "./output/表格3.1.1_发病数和发病率_逐年.xlsx")
yearly_incidence_county

if (!"发病率_区县_年" %in% names(wb1)) {addWorksheet(wb1, "发病率_区县_年")}
writeData(wb1, sheet = "发病率_区县_年", x = yearly_incidence_county)

# 读取地图
# nameCE_county$county_code  <-  as.character(nameCE_county$county_code)
map_data <- st_read("./data/origin/新疆维吾尔自治区.json") %>%
  st_transform(crs = 4326) %>%
  mutate(NAME = name, GB1999 = as.character(adcode), city_code = substr(GB1999,1,4)) 
# %>%
#   left_join(nameCE_county, by = c("adcode" = "county_code"))



# 检查哪些几何对象是无效的
invalid_index <- which(!st_is_valid(st_geometry(map_data)))

if (length(invalid_index) > 0) {
  cat("发现无效几何对象，尝试修复...\n")

  # 对无效几何对象进行修复
  map_data$geometry[invalid_index] <- st_make_valid(st_geometry(map_data)[invalid_index])

  # 再次检查是否仍然有无效几何对象
  invalid_after_first_fix <- which(!st_is_valid(st_geometry(map_data)))

  if (length(invalid_after_first_fix) > 0) {
    cat("仍有无效几何对象，尝试使用缓冲处理修复...\n")

    # 使用 st_buffer 修复剩余的无效几何对象
    buffer_dist <- 0.000001  # 使用一个非常小的正值避免重复边问题
    map_data$geometry[invalid_after_first_fix] <- st_buffer(
      st_geometry(map_data)[invalid_after_first_fix],
      dist = buffer_dist
    )

    # 再次检查是否仍然有无效几何对象
    invalid_after_buffer <- which(!st_is_valid(st_geometry(map_data)))

    if (length(invalid_after_buffer) > 0) {
      cat("仍有无效几何对象，尝试简化几何对象...\n")

      # 尝试简化几何对象
      simplify_tolerance <- 0.0001  # 根据需要调整简化容差
      map_data$geometry[invalid_after_buffer] <- st_simplify(
        st_geometry(map_data)[invalid_after_buffer],
        dTolerance = simplify_tolerance
      )

      # 最后再检查一次有效性
      final_invalid <- which(!st_is_valid(st_geometry(map_data)))
      if (length(final_invalid) > 0) {
        warning("仍有无效几何对象无法自动修复，请手动检查这些对象。\n")
      } else {
        cat("所有几何对象已成功修复。\n")
      }
    } else {
      cat("所有几何对象已通过缓冲处理成功修复。\n")
    }
  } else {
    cat("所有几何对象已通过 st_make_valid 成功修复。\n")
  }
} else {
  cat("所有几何对象都是有效的，无需修复。\n")
}

# 确保 map_data 的几何列是最新的
map_data <- st_as_sf(map_data)


# 检查哪些几何对象是无效的
invalid_index <- which(!st_is_valid(st_geometry(map_data)))


# 打印结果以确认
print(head(map_data))
map_data$adcode %>% unique() %in% "654223" %>% sum()
# class(map_data)



# 确保 county_code 和 GB1999 的类型一致
yearly_incidence_county <- yearly_incidence_county %>%
  mutate(county_code = as.character(county_code))

# 合并数据框
map_data <- map_data %>%
  filter(!city_code %in% c("6590", "6500"))

yearly_incidence_county$county_code %>% unique()
map_data$GB1999 %>% unique()
map_data$NAME %>% unique()

map_data$city_code %>% unique()

# 合并前检查county_code
# 提取唯一的 county_code 和 GB1999 值
county_codes <- unique(yearly_incidence_county$county_code)
gb1999_codes <- unique(map_data$GB1999)

# 找出相同元素（交集）
common_codes <- intersect(county_codes, gb1999_codes)

# 找出不同元素
# 在 yearly_incidence_county 中但不在 map_data 中的元素
only_in_county <- setdiff(county_codes, gb1999_codes)

# 在 map_data 中但不在 yearly_incidence_county 中的元素
only_in_gb1999 <- setdiff(gb1999_codes, county_codes)

# 输出结果
cat("相同的元素（交集）:\n")
print(common_codes)

cat("\n仅存在于 大疫情网 中的元素:\n")
print(only_in_county)

cat("\n仅存在于 地图中的元素:\n")
print(only_in_gb1999)

# 创建一个包含比较结果的数据框
# 确定最长的向量长度
max_length <- max(length(common_codes), length(only_in_county), length(only_in_gb1999))

# 将较短的向量扩展到最大长度，使用 NA 填充
common_codes_padded <- c(common_codes, rep(NA, max_length - length(common_codes)))
only_in_county_padded <- c(only_in_county, rep(NA, max_length - length(only_in_county)))
only_in_gb1999_padded <- c(only_in_gb1999, rep(NA, max_length - length(only_in_gb1999)))
# 【有个问题】"654223""654203"

# 创建一个包含比较结果的数据框
comparison_results <- data.frame(
  共有 = common_codes_padded,
  监测数据才有 = only_in_county_padded,
  地图才有 = only_in_gb1999_padded,
  stringsAsFactors = FALSE
)

# 将结果写入CSV文件
writexl::write_xlsx(comparison_results, "./output/表格3.1_地图和数据无法匹配的区县.xlsx")


yearly_incidence_county$city_code %>% unique()
map_data$city_code %>% unique()

aa <- yearly_incidence_county$county_code %>% unique()
bb <- map_data$adcode %>% unique()

common_codes <- intersect(aa, bb)

setdiff(aa, bb);
setdiff(bb, aa)

#合并
merged_data <- yearly_incidence_county %>%
    mutate(county_code = ifelse(county_code == "652403", "652423", county_code))%>%
  mutate(county_code = ifelse(county_code == "654203", "654223", county_code))%>%
  left_join(map_data, by = c("county_code" = "GB1999","city_code")) %>%
  dplyr::select(-name,-childrenNum,-adcode) %>%
  left_join(nameCE[c("nameC","nameE")], by = c("地区" = "nameC")) %>% 
  left_join(rate_year_county, by = c("地区","region", "year","county_code" = "地区编码")) %>% 
  rename(pop = value)

# merged_data_city <- yearly_incidence_city %>%
#   mutate(county_code = ifelse(county_code == "652403", "652423", county_code))%>%
#   mutate(county_code = ifelse(county_code == "654203", "654223", county_code))%>%
#   left_join(map_data, by = c("county_code" = "GB1999","city_code")) %>%
#   dplyr::select(-name,-childrenNum,-adcode) %>%
#   left_join(nameCE[c("nameC","nameE")], by = c("地区" = "nameC")) %>% 
#   left_join(rate_year_county, by = c("地区","region", "year","county_code" = "地区编码")) %>% 
#   rename(pop = value)
saveRDS(merged_data, "./data/processed/merged_data.rds")

names(merged_data)
head(merged_data)


# 3.2 绘图：每年发病数和发病率--------


# time_periods <- list(
#   "All Years" = range(data$year),
#   "2005-2009" = c(2005, 2009),
#   "2010-2015" = c(2010, 2015),
#   "2016-2019" = c(2016, 2019),
#   "2020-2023" = c(2020, 2023)
# )

# 创建一个空列表来存储每年的时间段
time_periods1 <- list()

# 使用for循环为每一年创建一个条目
for (year in 2005:2023) {
  time_periods1[[as.character(year)]] <- year}

time_periods <- time_periods1 #重要，分年份！！
# 重要：运行则是每年来；不运行则是分阶段

# 创建一个函数来计算每个时间段的平均发病率
calculate_incidence_rate <- function(data, years) {
  data %>%
    filter(year %in% years) %>%
    group_by(county_code, NAME, nameE) %>%
    # summarise(avg_incidence_rate = mean(incidence_rate, na.rm = TRUE), .groups = 'drop')
    summarise(
      total_cases = sum(count, na.rm = TRUE),      # 时段总发病数
      total_pop = sum(pop, na.rm = TRUE),          # 时段总人口数
      avg_incidence_rate = (total_cases / total_pop) * 100000,  # 每10万人发病率
      .groups = 'drop'
    ) %>%
    select(-total_cases, -total_pop)
}
# 
#     group_by(county_code, NAME) %>%
#     summarise(
#       total_cases = sum(count, na.rm = TRUE),      # 时段总发病数
#       total_pop = sum(pop, na.rm = TRUE),          # 时段总人口数
#       avg_incidence_rate = (total_cases / total_pop) * 100000,  # 每10万人发病率
#       .groups = 'drop'
#     ) %>%
#     select(-total_cases, -total_pop)

# 计算不同时间段的平均发病率
incidence_rates <- lapply(time_periods, calculate_incidence_rate, data = merged_data)
names(incidence_rates) <- names(time_periods)

# 将 incidence_rates 列表中的每个数据框合并到一个单独的数据框中，并且为每个数据框添加一个标识时间段的列
# 为每个子数据框添加时间段标识，并合并所有子数据框
combined_data <- map_df(incidence_rates, ~ .x %>% 
                          mutate(time_period = cur_group()), .id = "time_period")

# 将 time_period 列的值从数字转换回原始的时间段名称
combined_data <- combined_data %>%
  mutate(time_period = factor(time_period, levels = names(incidence_rates)))

# 查看合并后的数据框前几行
print(head(combined_data))
writexl::write_xlsx(combined_data, "./output/发病率_逐年.xlsx")

# wb1 <- createWorkbook()  # 新建空工作簿（无任何 Sheet）
if (!"发病率_逐年" %in% names(wb1)) {addWorksheet(wb1, "发病率_逐年")}
writeData(wb1, sheet = "发病率_逐年", x = combined_data)


# 3.3 绘图：分阶段发病数和发病率 --------

time_periods <- list(
  all_years = 2005:2023,
  y2005_2009 = 2005:2009,
  y2010_2015 = 2010:2015,
  y2016_2019 = 2016:2019,
  y2020_2023 = 2020:2023
)

# 创建一个函数来计算每个时间段的平均发病率
calculate_incidence_rate <- function(data, years) {
  data %>%
    filter(year %in% years) %>%
    group_by(county_code, NAME, nameE) %>%
    # summarise(avg_incidence_rate = mean(incidence_rate, na.rm = TRUE), .groups = 'drop')
    summarise(
      total_cases = sum(count, na.rm = TRUE),      # 时段总发病数
      total_pop = sum(pop, na.rm = TRUE),          # 时段总人口数
      avg_incidence_rate = (total_cases / total_pop) * 100000,  # 每10万人发病率
      .groups = 'drop'
    ) %>%
    select(-total_cases, -total_pop)  # 可选：移除中间计算列
}

# 计算不同时间段的平均发病率
incidence_rates <- lapply(time_periods, calculate_incidence_rate, data = merged_data)
names(incidence_rates) <- names(time_periods)

# 将 incidence_rates 列表中的每个数据框合并到一个单独的数据框中，并且为每个数据框添加一个标识时间段的列
combined_data <- map_df(incidence_rates, ~ .x %>%
                          mutate(time_period = cur_group()), .id = "time_period")

# 将 time_period 列的值从数字转换回原始的时间段名称
combined_data <- combined_data %>%
  mutate(time_period = factor(time_period, levels = names(incidence_rates)))
combined_data$time_period %>% unique()

# wb1 <- createWorkbook()  # 新建空工作簿（无任何 Sheet）
if (!"发病率_分阶段" %in% names(wb1)) {addWorksheet(wb1, "发病率_分阶段")}
writeData(wb1, sheet = "发病率_分阶段", x = combined_data)

combined_data$avg_incidence_rate %>% summary()
combined_data$county_code  %>% unique()

# 查看合并后的数据框前几行
print(head(combined_data))
summary(combined_data$avg_incidence_rate)

# 定义全局的 color_scale 范围
global_color_scale <- seq(0, 300, by = 50)

# 定义色板
# 定义颜色
display.brewer.all(n=10, exact.n=FALSE)

palette <- c(brewer.pal(9, "OrRd")[c(1:7)])


# 绘制地图的函数 
plot_incidence_map_tmap <- function(rate_data, title, global_color_scale, palette) {
  # 合并地图数据与发病率数据
  rate_map_data <- map_data %>%
    left_join(rate_data, by = c("GB1999" = "county_code"))
  # 创建 tmap 地图
  tm_shape(rate_map_data) +
    tm_polygons(col = "avg_incidence_rate",
                style = "fixed",
                palette = palette,
                breaks = global_color_scale,
                title = "Incidence Rate",
                textNA = "Data Missing",  # 图例中显示缺失值标签
                na.color = "white") +
    tm_layout(
      main.title = title,
      main.title.position = "center",
      legend.position = c("right", "bottom"),
      legend.text.size = 0.8,
      legend.title.size = 0.9,
      frame = FALSE
    ) +
    tm_compass(position = c("left", "top"), size = 2) +  # 添加指北针
    tm_scale_bar(position = c("left", "bottom"))        # 添加比例尺
}

# 确保目录存在，如果不存在则创建
dir.create("figs/逐年发病地图", showWarnings = FALSE, recursive = TRUE)

# 分别绘制不同时间段的地图并保存图片
plots <- lapply(names(incidence_rates), function(period) {
  p <- plot_incidence_map_tmap(
    incidence_rates[[period]],
    title = paste("Incidence Rate for", period),
    global_color_scale = global_color_scale,
    palette = palette
  )

  filename <- file.path("figs/逐年发病地图", paste0("incidence_rate_", period, ".png"))
  tmap_save(p, filename = filename, dpi = 300)})

# 3.4 空间自相关：分区县--------

# 3.4.1 空间自相关：2023年

# 过滤出2023年的数据
year_2023_data <- yearly_incidence_county %>%
  filter(year == 2023) %>%
  select(county_code, incidence_rate)

# 确保所有区县都在数据中，并用0填充缺失值
year_2023_map_data <- map_data %>%
  left_join(year_2023_data, by = c("GB1999" = "county_code")) %>%
  mutate(incidence_rate = coalesce(incidence_rate, 0L))

# 创建邻居关系
nb_q <- poly2nb(st_geometry(map_data), queen = TRUE)

# 创建权重列表，并明确指定 zero.policy = TRUE
lw_q <- nb2listw(nb_q, style = "W", zero.policy = TRUE)
length(lw_q$neighbours)

# 执行全局莫兰检验
global_moran_result <- moran.test(year_2023_map_data$incidence_rate, listw = lw_q,
                                  randomisation = FALSE, zero.policy = TRUE)

# 查看全局莫兰检验结果
print(global_moran_result)

# LISA结果
local_moran_result <- localmoran(year_2023_map_data$incidence_rate, listw = lw_q,
                                 zero.policy = TRUE)

# 整理局部莫兰检验结果为数据框
local_moran_df <- data.frame(
  county_code = year_2023_map_data$GB1999,
  name = year_2023_map_data$name,
  incidence_rate = year_2023_map_data$incidence_rate,
  Ii = local_moran_result[, 1],
  E_Ii = local_moran_result[, 2],
  Var_Ii = local_moran_result[, 3],
  Z_Ii = local_moran_result[, 4],
  Pr_Ii = local_moran_result[, 5]
)

# 查看局部莫兰检验结果的前几行
print(head(local_moran_df))

# 3.4.2 空间自相关：2005-2023年

# 步骤 1：定义函数和初始化

# 创建邻居关系
nb_q <- poly2nb(st_geometry(map_data), queen = TRUE)

# 创建权重列表，并明确指定 zero.policy = TRUE
lw_q <- nb2listw(nb_q, style = "W", zero.policy = TRUE)

# 定义一个函数来执行全局莫兰检验
perform_global_moran <- function(data) {
  moran.test(data$incidence_rate, listw = lw_q, randomisation = FALSE, zero.policy = TRUE)
}

# 定义一个函数来执行局部莫兰检验（LISA），并添加聚集类型分类
perform_local_moran <- function(data) {
  local_moran_result <- localmoran(data$incidence_rate, listw = lw_q, zero.policy = TRUE)
  
  # 创建数据框
  moran_df <- data.frame(
    county_code = data$GB1999,
    name = data$name,
    incidence_rate = data$incidence_rate,
    Ii = local_moran_result[, 1],
    E_Ii = local_moran_result[, 2],
    Var_Ii = local_moran_result[, 3],
    Z_Ii = local_moran_result[, 4],
    Pr_Ii = local_moran_result[, 5]
  )
  
  # 添加聚集类型分类
  mean_rate <- mean(data$incidence_rate, na.rm = TRUE)
  moran_df$cluster_type <- ifelse(moran_df$incidence_rate > mean_rate,
                                  ifelse(moran_df$Ii > 0, "High-High", "High-Low"),
                                  ifelse(moran_df$Ii > 0, "Low-High", "Low-Low"))
  
  moran_df}

# 步骤 2：循环处理每年的数据

# 获取所有年份
years <- unique(yearly_incidence_county$year)

# 初始化结果列表
global_moran_results <- list()
local_moran_results <- list()

# 循环每年的数据并进行检验
for (yr in years) {
  # 过滤出当年的数据，并确保所有区县都在数据中，并用0填充缺失值
  yearly_data <- yearly_incidence_county %>%
    filter(year == yr) %>%
    select(county_code, incidence_rate)
  
  yearly_map_data <- map_data %>%
    left_join(yearly_data, by = c("GB1999" = "county_code")) %>%
    mutate(incidence_rate = coalesce(incidence_rate, 0L))
  
  # 执行全局莫兰检验
  global_moran_results[[as.character(yr)]] <- perform_global_moran(yearly_map_data)
  
  # 执行局部莫兰检验（LISA）
  local_moran_results[[as.character(yr)]] <- perform_local_moran(yearly_map_data)
}

# 步骤 3：整理全局莫兰检验结果
# 将全局莫兰检验结果整理成数据框
global_moran_df <- data.frame(
  year = as.numeric(names(global_moran_results)),
  moran_I = sapply(global_moran_results, function(x) x$estimate["Moran I statistic"]),
  p_value = sapply(global_moran_results, function(x) x$p.value),
  stringsAsFactors = FALSE
)

# 查看全局莫兰检验结果
print(head(global_moran_df))

wb2 <- createWorkbook()  # 新建空工作簿（无任何 Sheet）
if (!"全局莫兰_逐年" %in% names(wb2)) {addWorksheet(wb2, "全局莫兰_逐年")}
writeData(wb2, sheet = "全局莫兰_逐年", x = global_moran_df)

# 步骤 4：整理局部莫兰检验结果

# 初始化结果列表
cluster_summary <- list()

# 整理局部莫兰检验结果
for (yr in names(local_moran_results)) {
  yearly_local_moran <- local_moran_results[[yr]] %>% 
    filter(Pr_Ii < 0.05)
  
  cluster_counts <- yearly_local_moran %>%
    group_by(cluster_type) %>%
    summarise(count = n(), 
              counties = paste(name, collapse = ", ")) %>%
    ungroup() %>%
    complete(cluster_type = c("High-High", "Low-Low", "High-Low", "Low-High"), fill = list(count = 0, counties = ""))
  
  cluster_summary[[yr]] <- cluster_counts
}
str(local_moran_results)

# 3.4.2.1 合并LISA的所有指标的结果
# 初始化一个空列表来保存调整后的数据框
adjusted_dataframes <- list()

# 遍历 local_moran_results 列表，为每个数据框添加 'year' 列，并存入 adjusted_dataframes 列表
for (year in names(local_moran_results)) {
  df <- local_moran_results[[year]]
  df$year <- as.numeric(year) # 添加 'year' 列
  adjusted_dataframes[[year]] <- df
}

# 使用 bind_rows() 将所有调整后的数据框合并成一个大的数据框
combined_df <- bind_rows(adjusted_dataframes) #%>% filter(Pr_Ii < 0.05);2.8要作图；不能这么设定
writexl::write_xlsx(combined_df, "./output/表格3.4.2_LISA_历年.xlsx")

head(combined_df)
saveRDS(combined_df, "./output/LISA结果_历年.rds")


if (!"LISA_逐年" %in% names(wb2)) {addWorksheet(wb2, "LISA_逐年")}
writeData(wb2, sheet = "LISA_逐年", x = combined_df)


# # 3.4.2.2 合并LISA的空间聚集分析结果
# 合并结果
cluster_summary_df <- bind_rows(cluster_summary, .id = "year")
saveRDS(cluster_summary_df, "./output/LISA结果_计数_历年.rds")


# 调整列顺序以便查看
cluster_summary_df <- cluster_summary_df %>%
  select(year, cluster_type, count, counties)

# 查看局部莫兰检验结果
print(head(cluster_summary_df))

if (!"LISA_逐年_计数" %in% names(wb2)) {addWorksheet(wb2, "LISA_逐年_计数")}
writeData(wb2, sheet = "LISA_逐年_计数", x = cluster_summary_df)


# 步骤 5：保存结果
global_moran_df$year %>% unique()
cluster_summary_df$year %>% unique()

writexl::write_xlsx(global_moran_df, "./output/表格3.4.1_全局莫兰.xlsx")
writexl::write_xlsx(cluster_summary_df, "./output/表格3.4.2_局部莫兰_LISA_计数.xlsx")

# 步骤 6：将局部莫兰指数转宽数据
cluster_summary_df1 <- cluster_summary_df %>% dplyr::select(-counties)
cluster_summary_df1
cluster_summary_df2 <- cluster_summary_df1 %>% 
  pivot_wider(names_from = year, values_from = count) 
writexl::write_xlsx(cluster_summary_df2, "./output/表格3.4.3_LISA_宽数据_历年.xlsx")

if (!"LISA_逐年_计数_宽数据" %in% names(wb2)) {addWorksheet(wb2, "LISA_逐年_计数_宽数据")}
writeData(wb2, sheet = "LISA_逐年_计数_宽数据", x = cluster_summary_df2)


# 步骤 7：将历年聚集类型绘图

cluster_summary_df1$year <- as.numeric(cluster_summary_df1$year)

# 折线图：不同类型聚集区县数量
ggplot(cluster_summary_df1, aes(x = year, y = count, color = cluster_type, group = cluster_type)) +
  geom_line(size = 1) + # 添加线条
  geom_point(size = 3) + # 添加点，使每个数据点更明显
  labs(
    title = "Local Moran's I Analysis of Brucellosis Incidence in Xinjiang (2005-2023)",
    x = "Year",
    y = "Number of Counties",
    color = "Cluster Type"
  ) +
  theme_minimal() + # 使用简洁的主题
  scale_color_manual(values = RColorBrewer::brewer.pal(4, "Set3")) + # 使用 Set3 色系
  theme(
    plot.title = element_text(hjust = 0.5), # 居中标题
    axis.text.x = element_text(angle = 45, hjust = 1) # 调整x轴标签角度以便阅读
  )

ggsave("./figs/图3.4.2 折线图_LISA_历年.png")

# 3.4.3 空间自相关：分4阶段--------
head(map_data)

time_periods <- list(
  "All Years" = c(2005, 2023),
  "2005-2009" = c(2005, 2009),
  "2010-2015" = c(2010, 2015),
  "2016-2019" = c(2016, 2019),
  "2020-2023" = c(2020, 2023)
)

# 创建邻居关系和权重列表
nb_q <- poly2nb(st_geometry(map_data), queen = TRUE)
lw_q <- nb2listw(nb_q, style = "W", zero.policy = TRUE)
length(lw_q$neighbours)

# 初始化结果列表
global_moran_results <- list()
local_moran_results <- list()

# 循环每个时间段进行全局和局部莫兰检验
for (period_name in names(time_periods)) {
  years_in_period <- time_periods[[period_name]]
  
  # 聚合每年的数据，并计算平均发病率
  period_map_data <- merged_data %>%
    filter(year %in% years_in_period) %>%
    group_by(county_code, NAME,nameE) %>%
    # summarise(incidence_rate = mean(incidence_rate, na.rm = TRUE), .groups = 'drop')
  # filter(year %in% years) %>%
  #   group_by(county_code, NAME,nameE) %>%
    summarise(
      total_cases = sum(count, na.rm = TRUE),      # 时段总发病数
      total_pop = sum(pop, na.rm = TRUE),          # 时段总人口数
      incidence_rate = (total_cases / total_pop) * 100000,  # 每10万人发病率
      .groups = 'drop'
    ) %>%
    mutate(incidence_rate = ifelse(is.nan(incidence_rate), 0, incidence_rate)) %>%
    select(-total_cases, -total_pop)  # 可选：移除中间计算列

  # merged_data <- yearly_incidence_county %>%
  #   mutate(county_code = ifelse(county_code == "652403", "652423", county_code))%>%
  # mutate(county_code = ifelse(county_code == "654203", "654223", county_code))%>%
  # left_join(map_data, by = c("county_code" = "GB1999","city_code")) %>%
  # dplyr::select(-name,-childrenNum,-adcode) %>%
  # left_join(nameCE[c("nameC","nameE")], by = c("地区" = "nameC")) %>% 
  # left_join(rate_year_county, by = c("地区","region", "year","county_code" = "地区编码")) %>% 
  # rename(pop = value)
  
  # 确保所有区县都在数据中，并用0填充缺失值
  # period_map_data <- map_data %>%
  #   left_join(period_data, by = c("GB1999" = "county_code")) %>%
  #   mutate(incidence_rate = coalesce(incidence_rate, 0L))
  
  # 执行全局莫兰检验
  global_moran_results[[period_name]] <- moran.test(period_map_data$incidence_rate, lw_q, zero.policy = TRUE)
  
  # 执行局部莫兰检验（LISA），并添加聚集类型分类
  local_moran_result <- localmoran(period_map_data$incidence_rate, lw_q, zero.policy = TRUE)
  
  # 将局部莫兰检验结果与原始数据结合，并添加聚集类型分类
  local_moran_df <- data.frame(
    period = period_name,
    name = period_map_data$nameE,
    GB1999 = period_map_data$county_code,
    incidence_rate = period_map_data$incidence_rate,
    Ii = local_moran_result[, 1],
    E_Ii = local_moran_result[, 2],
    Var_Ii = local_moran_result[, 3],
    Z_Ii = local_moran_result[, 4],
    Pr_Ii = local_moran_result[, 5]
  ) %>%
    mutate(cluster_type = case_when(
      incidence_rate > mean(incidence_rate, na.rm = TRUE) & Ii > 0 ~ "High-High",
      incidence_rate > mean(incidence_rate, na.rm = TRUE) & Ii <= 0 ~ "High-Low",
      incidence_rate <= mean(incidence_rate, na.rm = TRUE) & Ii > 0 ~ "Low-High",
      incidence_rate <= mean(incidence_rate, na.rm = TRUE) & Ii <= 0 ~ "Low-Low"
    ))
  
  local_moran_results[[period_name]] <- local_moran_df
}

local_moran_results1 <- map_dfr(names(local_moran_results), ~ {
  period_local_moran1 <- local_moran_results[[.x]]})
local_moran_results1 %>% filter(Pr_Ii < 0.05)

if (!"LISA_分阶段" %in% names(wb2)) {addWorksheet(wb2, "LISA_分阶段")}
writeData(wb2, sheet = "LISA_分阶段", x = local_moran_results1)

# 全局莫兰
global_moran_df <- data.frame(
  period = names(global_moran_results),
  moran_I = sapply(global_moran_results, function(x) x$estimate["Moran I statistic"]),
  p_value = sapply(global_moran_results, function(x) x$p.value),
  stringsAsFactors = FALSE
)

if (!"全局莫兰_分阶段" %in% names(wb2)) {addWorksheet(wb2, "全局莫兰_分阶段")}
writeData(wb2, sheet = "全局莫兰_分阶段", x = global_moran_df)

# LISA结果

cluster_summary_df <- map_dfr(names(local_moran_results), ~ {
  period_local_moran <- local_moran_results[[.x]] %>% 
    filter(Pr_Ii < 0.05)
  
  cluster_counts <- period_local_moran %>%
    group_by(cluster_type) %>%
    summarise(count = n(), 
              counties = paste(name, collapse = ", "), .groups = 'drop') %>%
    complete(cluster_type = c("High-High", "Low-Low", "High-Low", "Low-High"), 
             fill = list(count = 0, counties = ""))
  
  cluster_counts$period <- .x
  
  cluster_counts
})

# 查看结果
print(global_moran_df)
print(head(cluster_summary_df))

# 保存结果
writexl::write_xlsx(global_moran_df, "./output/表格3.4.3_空间自相关_全局_分阶段_.xlsx")
writexl::write_xlsx(cluster_summary_df, "./output/表格3.4.3_空间自相关_局部_分阶段.xlsx")

if (!"LISA_分阶段_计数" %in% names(wb2)) {addWorksheet(wb2, "LISA_分阶段_计数")}
writeData(wb2, sheet = "LISA_分阶段_计数", x = cluster_summary_df)

# 3.4.4 分阶段LISA地图
library(tmap)
library(dplyr)

# 创建“LISA地图”文件夹（如果不存在）
# lisa_map_dir <- "./figs/LISA地图"
lisa_map_dir <- "./主图/LISA地图"
if (!dir.exists(lisa_map_dir)) {
  dir.create(lisa_map_dir, showWarnings = FALSE, recursive = TRUE)
}

# # 循环每个时间段进行绘图
# for (period_name in names(local_moran_results)) {
#   period_local_moran <- local_moran_results[[period_name]]
# 
#   # 将局部莫兰检验结果与原始地图数据结合
#   map_data_lisa <- map_data %>%
#     left_join(period_local_moran, by = c("GB1999" = "GB1999"))
# 
#   # 创建 tmap 图层
#   tm <- tm_shape(map_data_lisa) +
#     tm_polygons(col = "cluster_type",
#                 title = "Cluster Type",
#                 palette = c("High-High" = "red",
#                             "Low-Low" = "blue",
#                             "High-Low" = "lightblue",
#                             "Low-High" = "pink"),
#                 border.col = "grey50") +
#     tm_layout(title = paste("Local Moran's I for", period_name),
#               legend.position = c("right", "bottom"))
# 
#   # 保存地图到文件
#   filename <- file.path(lisa_map_dir, paste0("分阶段空间自相关_局部莫兰检验_", period_name, ".png"))
#   tmap_save(tm, filename = filename, width = 12, height = 8, dpi = 300)
# }

# 
# # # 定义颜色映射（使用 Set3 色系）
# # if (TRUE) {  colors <- c("#8DD3C7", "#FFFFB3", "#BEBADA", "#FB8072", "#80B1D3", "#FDB462", "#B3DE69", "#D9D9D9", "#BC80BD", "#CCEBC5", "#FFED6F")
# # print(colors)}
# # palette <- "Reds"
# palette <- c("#80B1D3","#FFED6F","#FDB462",  "#BEBADA","#FB8072", "#D9D9D9" )
# brewer.pal(9, "Set3")[c(1,4)]

# cluster_colors <- RColorBrewer::brewer.pal(4, "Set3")
# cluster_colors <- c("#FF0000", "#0000FF",  "#00FFFF", "#FFA500", "#D3D3D3")
# cluster_colors <- list(
#   # 方案1：传统红蓝对比（学术论文常用）
#   cluster_colors1 = c(HH = "#FF0000", LL = "#0000FF", HL = "#00FFFF", LH = "#FFA500", NS = "#D3D3D3"),
# 
#   # 方案2：色盲友好配色（基于ColorBrewer）
#   cluster_colors2 = c(HH = "#E66101", LL = "#5E3C99", HL = "#B2ABD2", LH = "#FDB863", NS = "#F0F0F0"),
# 
#   # 方案3：自然渐变风格
#   cluster_colors3 = c(HH = "#D53E4F", LL = "#3288BD", HL = "#66C2A5", LH = "#FDAE61", NS = "#EDEDED"),
# 
#   # 方案4：高对比度组合
#   cluster_colors4 = c(HH = "#E41A1C", LL = "#377EB8", HL = "#984EA3", LH = "#FF7F00", NS = "#CCCCCC"),
# 
#   # 方案5：柔和色调（适合演示场景）
#   cluster_colors5 = c(HH = "#FF6B6B", LL = "#4ECDC4", HL = "#8395A7", LH = "#FF9F43", NS = "#F5F5F5")
# )
# cluster_colors <-  c("#FB8072","#8DD3C7","#FDB462","#BEBADA" )

# cluster_colors <-  c("#FF0000","#0000FF","#00FFFF","#FFA500" )
# names(cluster_colors) <- c("High-High", "Low-Low", "High-Low", "Low-High")

# # 循环每个时间段进行绘图
# for (period_name in names(local_moran_results)) {
#   period_local_moran <- local_moran_results[[period_name]]
#   
#   # 将局部莫兰检验结果与原始地图数据结合
#   map_data_lisa <- map_data %>%
#     left_join(period_local_moran, by = c("GB1999" = "GB1999"))
#   
#   # 创建 tmap 图层
#   tm <- tm_shape(map_data_lisa) +
#     tm_polygons(col = "cluster_type", 
#                 title = "Cluster Type",
#                 # palette = cluster_colors,  # 使用 Set3 色系
#                 border.col = "grey50") +
#     tm_layout(title = paste("Local Moran's I for", period_name),
#               legend.position = c("right", "bottom"),
#               main.title = "Local Moran's I Analysis",
#               main.title.position = "center")
#   
#   # 保存地图到文件
#   filename <- file.path(lisa_map_dir, paste0("local_moran_", period_name, ".png"))
#   tmap_save(tm, filename = filename, width = 12, height = 8, dpi = 300)
# }
# names(local_moran_results)
# local_moran_results = "2005-2009"
# period_name = "2005-2009"

# colors <- c("red",  "#FFB3BA", "#B3C6FF","blue")
# tm_polygons(..., palette = colors)
# 
# for (period_name in names(local_moran_results)) {
#   period_local_moran <- local_moran_results[[period_name]] %>% 
#     dplyr::select(GB1999, Pr_Ii, cluster_type)
#   
#   # 将局部莫兰检验结果与原始地图数据结合
#   map_data_lisa <- map_data %>%
#     left_join(period_local_moran, by = c("GB1999" = "GB1999")) 
#   
#   # dev.new()
#   # 创建 tmap 图层
#   tm <- tm_shape(map_data_lisa) +
#     tm_polygons(
#       col = "cluster_type",  ## map_data_lisa$cluster_type %>% unique()
#       palette = colors,
#       title = "Cluster Type",
#       # textNA = "Data Missing",  # 图例中显示缺失值标签
#       na.color = "white",  # 缺失值区域颜色设为白色
#       border.col = "grey50"
#     ) +
#     # 添加指北针（右上角，箭头样式）
#     tm_compass(
#       type = "arrow", 
#       position = c("right", "top"),
#       size = 2,  # 控制尺寸
#       text.size = 0.8  # 文字大小
#     ) +
#     # 添加比例尺（左下角，公制单位）
#     tm_scale_bar(
#       # position = c("left", "bottom"), #breaks和width不能一起使用，width被忽略
#       width = 0.2,  # 比例尺长度占地图宽度的比例
#       text.size = 0.7,  # 文字大小
#       # breaks = c(0, 100, 200),  # 主刻度值（单位：公里）
#       text.color = "black"
#     ) +
#     tm_layout(
#       # title = paste("Local Moran's I for", period_name),
#       legend.position = c("right", "bottom"),
#       # main.title = "Local Moran's I Analysis",
#       main.title = paste("Local Moran's I for", period_name),
#       main.title.position = "center"
#     )
#   
#   # 保存地图到文件
#   filename <- file.path(lisa_map_dir, paste0("local_moran_", period_name, ".png"))
#   tmap_save(tm, filename = filename, dpi = 300)
# }

# 初始化空列表
map_list <- list()

# 需求一：显示完整类别名称的图例
for (period_name in names(local_moran_results)) {
  period_local_moran <- local_moran_results[[period_name]] %>% 
    dplyr::select(GB1999, Pr_Ii, cluster_type,name)
  
  # 数据处理：创建全称因子
  map_data_lisa <- map_data %>%
    left_join(period_local_moran, by = "GB1999") %>% 
    mutate(
      cluster_category = case_when(
        is.na(Pr_Ii) | Pr_Ii >= 0.05 ~ "Not Significant",
        TRUE ~ as.character(cluster_type)
      ),
      cluster_category = factor(
        cluster_category,
        levels = c("High-High",  "High-Low", "Low-High", "Low-Low","Not Significant"),#"High-Low", "Low-High", 
        labels = c("High-High",   "High-Low", "Low-High","Low-Low","Not Significant") #"High-Low", "Low-High",
      )
    )
  
  # 完整颜色映射
  full_colors <- c(
    "High-High" = "red",
    "High-Low" = "#FFB3BA",
    "Low-High" = "#B3C6FF",
    "Low-Low" = "blue",
    "Not Significant" = "white")
  
  # 创建地图
  tm <- tm_shape(map_data_lisa) +
    tm_layout(bg.color = "transparent") +
    tm_polygons(
      col = "cluster_category",
      palette = full_colors,
      title = "Cluster Type",
      border.col = "grey50",
      colorNA = NULL,
      textNA = NULL,
      legend.show = F
    ) +
    tm_layout(
      main.title = paste( period_name),
      # legend.outside = FALSE
      main.title.position = "center",
      main.title.fontface = "bold",
      fontfamily = "Arial",
      # legend.position = c("right", "bottom"),
      # legend.bg.color = "white",
      # legend.bg.alpha = 0.8,
      # legend.frame = FALSE,
      # legend.text.size = 0.85,
      # legend.title.size = 1.0,
      legend.show = F,
      frame = FALSE,
      outer.margins = c(0.02, 0, 0.02, 0)
    )
  # 动态命名列表元素
  map_list[[period_name]] <- tm  # 用时间段名称作为列表键
  
  # 保存输出
  filename <- file.path(lisa_map_dir, paste0("full_legend_", period_name, ".png"))
  tmap_save(tm, filename = filename, dpi = 600)
}



map_list[[1]]
# 
# # 循环保存每张地图
# for(i in seq_along(map_list)) {
#   # 生成文件名（如："分阶段_LISA图1.png"）
#   filename <- sprintf("分阶段_LISA图%d.png", i) 
#   filepath <- file.path("./主图/LISA地图", filename)
#   
#   # 保存地图（适配ggplot2/tmap两种格式）
#   if(inherits(map_list[[i]], "ggplot")){
#     ggsave(filepath, plot = map_list[[i]], dpi = 600) # width = 10, height = 8,
#   } else if(inherits(map_list[[i]], "tmap")){
#     tmap::tmap_save(map_list[[i]], filename = filepath) #, width = 2400, height = 1800
#   }
# }


# 创建第20张图：比例尺、指北针和图例

legend_map <- tm_shape(map_data) +
  tm_polygons() +
  tm_layout(bg.color = "transparent") +
  tm_compass(#position = c("left", "top"),
    position = c(0.2,0.8),
    size = 1.5) +
  tm_scale_bar(
    # position = c("right", "top"),
    position = c(0.2,0.65),
    breaks = c(0,400,800),
    text.size = 0.7) +
  tm_add_legend(
    type = "fill",
    labels = c("High-High",  "High-Low", "Low-High","Low-Low", "Not Significant"),
    col = c("#d73027", "#fc8d59",  "#91bfdb","blue","white"),
    title = "Cluster Category",
    border.lwd = 0.5
  ) +
  tm_layout(
    legend.only = TRUE,
    legend.position = c(0.2, 0.2),
    legend.width = 1.5,
    legend.title.size = 1.2,
    legend.text.size = 0.8
  )

# 获取当前列表名称并重新排序
new_order <- c(setdiff(names(map_list), "All Years"), "All Years")

# 按新顺序重组列表
map_list <- map_list[new_order]

# 将 yearly_maps 和 legend_map 合并为一个列表
combined_list <- c(map_list,list(legend_map))

# 组合输出 (5列4行布局)
final_map <- tmap_arrange(
  combined_list,
  ncol = 3,  # 每行 5 张图
  nrow = 2,  # 每列 4 张图
  outer.margins = c(0.02, 0, 0.02, 0),
  asp = 0.6
)
final_map

tmap_save(
  final_map,
  filename = "./主图/LISA地图/分阶段_LISA图.png",
  # width = 45,     # 宽度 45cm（期刊通常接受 10-50cm）
  # height = 35,    # 高度 35cm
  # units = "cm",   # 单位厘米（期刊推荐）
  dpi = 600       # 分辨率 600 DPI（期刊要求 ≥ 300 DPI）
)


# # 需求二：仅显示High-High类别
# for (period_name in names(local_moran_results)) {
#   period_local_moran <- local_moran_results[[period_name]] %>% 
#     dplyr::select(GB1999, Pr_Ii, cluster_type)
#   
#   # 数据处理 
#   map_data_lisa <- map_data %>%
#     left_join(period_local_moran, by = "GB1999") %>% 
#     mutate(
#       cluster_category = if_else(
#         !is.na(Pr_Ii) & Pr_Ii < 0.05 & cluster_type == "High-High",
#         "High-High", 
#         NA_character_
#       ),
#       cluster_category = factor(cluster_category, levels = "High-High")
#     ) %>% mutate(adcode = as.character(adcode)) %>% 
#     left_join(name_county, by = c("adcode" = "county_code")) %>% 
#     mutate(name = nameE)
#   
#   # 创建标注层（仅包含显著High-High区域）
#   label_data <- map_data_lisa %>% 
#     filter(cluster_category == "High-High" & !is.na(name))
#   
#   # 可视化 
#   tm <- tm_shape(map_data_lisa) +
#     # 基础多边形
#     tm_polygons(
#       col = "cluster_category",
#       # palette = c("High-High" = "red"),
#       palette = c("High-High" = "#E41A1C"), # 使用更醒目的红色
#       title = "Hot Spots",
#       border.col = "grey50",
#       # colorNA = "white",
#       colorNA = NULL,  # 禁用默认缺失颜色
#       legend.show = TRUE
#     ) +
#     # 添加名称标注
#     tm_shape(label_data) +  # 使用过滤后的数据层
#     tm_text(
#       text = "name",
#       col = "black",
#       size = 1.3,        # 主标签大小 (期刊推荐相对大小)
#       fontface = "bold", # 加粗字体
#       alpha = 0.9,
#       # bg.color = "white",# 白色背景衬底
#       bg.alpha = 0.4,    # 背景透明度
#       auto.placement = 0.6,  # 增加自动布局间距
#       remove.overlap = TRUE,
#       just = "center"
#     ) +
#     # 其他地图元素
#     tm_compass(
#       type = "arrow", 
#       position = c("right", "top"),
#       # size = 2,
#       size = 3,
#       text.size = 1.0,   # 比例尺文字放大
#       # text.size = 0.8
#     ) +
#     tm_scale_bar(
#       width = 0.2,
#       text.size = 0.7,
#       text.color = "black"
#     ) +
#     tm_layout(
#       main.title = paste("High-High Clusters in", period_name),
#       # main.title = paste("High-High Clusters:", period_name),
#       main.title.size = 1.4,     # 主标题放大
#       main.title.position = "center",
#       fontfamily = "Arial",      # 使用期刊推荐字体
#       legend.title.size = 1.2,
#       legend.text.size = 1.0,
#       frame = FALSE
#     )
#   
#   # 保存输出
#   filename <- file.path(lisa_map_dir, paste0("highhigh_label_", period_name, ".png"))
#   tmap_save(tm, filename = filename, dpi = 300, width = 10, height = 8)
# }
# names(map_data_lisa)
# map_data_lisa$adcode %>% unique()
# map_data_lisa$cluster_type %>% unique()
# 
# 
# # 确认图片已保存到指定文件夹
# cat("LISA 地图已成功保存至：", lisa_map_dir, "\n")

# # 4、地区分布：分地市------
# 
# # 4.1. 2005-2023年每年的发病数和发病率（作图用）------
# data$county_code <- data$city_code
# 
# yearly_cases_county <- data %>%
#   filter(year >= 2004, year <= 2023) %>%
#   count(year, county_code, name = "cases")
# 
# yearly_incidence_county <- rate2 %>%
#   filter(region == "city") %>% #county
#   mutate(county_code = substr(county_code,1,4)) %>%
#   left_join(yearly_cases_county, by = c("year", "county_code")) %>%
#   select(-地区编码)
# writexl::write_xlsx(yearly_incidence_county, "./output/表格3.1.1_yearly_incidence_city.xlsx")
# 
# 
# # 读取地图
# map_data <- st_read("./data/origin/新疆维吾尔自治区.json")  %>%
#   mutate(GB1999 = as.character(GB1999))
# 
# # 确保 county_code 和 GB1999 的类型一致
# yearly_incidence_county <- yearly_incidence_county %>%
#   mutate(county_code = as.character(county_code))
# 
# # map_data <- map_data %>%
# #   mutate(GB1999 = as.character(GB1999)) %>%
# #   mutate(GB1999 = substr(GB1999,1,4))
# 
# # 合并数据框
# merged_data <- yearly_incidence_county #%>%
#   # left_join(map_data, by = c("county_code" = "GB1999"))
# 
# # 4.2 表格：分地州，不同时间段发病数和发病率--------
# time_periods <- list(
#   all_years = 2004:2023,
#   y2005_2015 = 2005:2015,
#   y2016_2020 = 2016:2020,
#   y2021_2023 = 2021:2023
# )
# 
# # 创建一个函数来计算每个时间段的平均发病率
# calculate_incidence_rate <- function(data, years) {
#   data %>%
#     filter(year %in% years) %>%
#     group_by(county_code) %>%#, NAME
#     summarise(avg_incidence_rate = mean(incidence_rate, na.rm = TRUE), .groups = 'drop')
# }
# 
# # 计算不同时间段的平均发病率
# incidence_rates <- lapply(time_periods, calculate_incidence_rate, data = merged_data)
# names(incidence_rates) <- names(time_periods)
# 
# # 将 incidence_rates 列表中的每个数据框合并到一个单独的数据框中，并且为每个数据框添加一个标识时间段的列
# # 为每个子数据框添加时间段标识，并合并所有子数据框
# combined_data <- map_df(incidence_rates, ~ .x %>% mutate(time_period = cur_group()), .id = "time_period")
# 
# # 将 time_period 列的值从数字转换回原始的时间段名称
# combined_data <- combined_data %>%
#   mutate(time_period = factor(time_period, levels = names(incidence_rates)))
# 
# # 查看合并后的数据框前几行
# print(head(combined_data))
# 
# writexl::write_xlsx(combined_data, "./output/表格4.2_incidence_rates_city_period.xlsx")
# 
# # 4.3 不绘制地州的地图了


# 5 新发感染区县增加的情况------

# 5.1 2005-2023年历年新发感染区县------
unique(yearly_incidence_county$city_code)
map_data$city_code %>% unique()


# 获取所有年份和区县的组合，并初始化状态列
full_grid <- expand.grid(year = seq(min(yearly_incidence_county$year), 
                                    max(yearly_incidence_county$year)),
                         county_code = unique(yearly_incidence_county$county_code)) %>%
  left_join(yearly_incidence_county, by = c("year", "county_code")) %>%
  mutate(status = "无病例报道")

# 更新状态列
update_status <- function(df) {
  df %>%
    group_by(county_code) %>%
    arrange(year) %>%
    mutate(
      cumulative_infections = cumsum(!is.na(cases) & cases > 0),
      status = case_when(
        !is.na(cases) & cases > 0 & cumulative_infections == 1 ~ "当年首次感染",
        !is.na(cases) & cases > 0 ~ "曾经感染过",
        is.na(cases) | cases == 0 ~ "无病例报道"
      )
    ) %>%
    ungroup() %>%
    mutate(status = ifelse(is.na(cumulative_infections), "从来没有感染过", status)) %>%
    select(-cumulative_infections)
}

full_grid <- update_status(full_grid)
writexl::write_xlsx(full_grid, "./output/表格5_emerging_infection_status.xlsx")



# 5.3.1 绘逐年地图 
# 绘逐年地图 
# plot_status_map <- function(rate_data, title, map_data) {
#   rate_map_data <- map_data %>%
#     left_join(rate_data, by = c("GB1999" = "county_code"))
# 
#   centroids <- st_centroid(st_geometry(rate_map_data))
#   rate_map_data$centroid_x <- st_coordinates(centroids)[,1]
#   rate_map_data$centroid_y <- st_coordinates(centroids)[,2]
# 
#   ggplot() +
#     geom_sf(data = rate_map_data, aes(fill = status)) +
#     scale_fill_manual(values = c("从来没有感染过" = "lightgreen",
#                                  "曾经感染过" = "orange",
#                                  "当年首次感染" = "red",
#                                  "无病例报道" = "white"),
#                       na.value = "gray") +
#     labs(title = title, fill = "Status") +
#     theme_void() +  # 使用 theme_void() 移除所有坐标轴和背景元素
#     theme(
#       plot.background = element_rect(fill = "white"),  # 设置背景为白色
#       plot.title = element_text(hjust = 0.5),          # 标题居中
#       legend.position = "bottom",                      # 将图例位置调整到下方
#       legend.box.background = element_rect(fill = "white", color = NA)  # 图例背景为白色
#     ) +
#     annotation_north_arrow(location = "tl", which_north = "true",
#                            style = north_arrow_fancy_orienteering) +
#     annotation_scale(location = "bl") +
#     geom_text(data = rate_map_data, aes(x = centroid_x, y = centroid_y, label = NAME), size = 3, vjust = -0.5)
# }
# 
# # 确保目录存在，如果不存在则创建
# dir.create("figs/逐年状态地图", showWarnings = FALSE, recursive = TRUE)
# 
# # 分别绘制不同年份的地图并保存图片
# plots <- lapply(unique(full_grid$year), function(yr) {
#   year_data <- full_grid %>% filter(year == yr)
# 
#   p <- plot_status_map(year_data, title = paste("Infection Status for", yr), map_data = map_data)
# 
#   filename <- file.path("figs/逐年状态地图", paste0("infection_status_", yr, ".png"))
#   ##ggsave(filename, plot = p, width = 10, height = 8, dpi = 300)
# 
#   p
# })
plot_status_map <- function(rate_data, title, map_data) {
  # 定义完整的状态水平（含NA处理）
  status_levels <- c("Never Infected", "Historically Infected", 
                     "Newly Infected", "No Case Reported", "Data Missing")
  
  # 数据预处理管道
  rate_map_data <- map_data %>%
    left_join(rate_data, by = c("GB1999" = "county_code")) %>%
    mutate(
      # 将中文标签转换为英文
      status = case_when(
        status == "从来没有感染过" ~ "Never Infected",
        status == "曾经感染过" ~ "Historically Infected",
        status == "当年首次感染" ~ "Newly Infected",
        status == "无病例报道" ~ "No Case Reported",
        TRUE ~ NA_character_
      ),
      # 处理合并产生的NA值
      status = ifelse(is.na(status), "Data Missing", status),
      # 转换为有序因子
      status = factor(status, levels = status_levels)
    )
  
  # 设置色板（使用Set3色系+补充灰色）
  color_palette <- c(RColorBrewer::brewer.pal(12, "Set3")[c(3,7,9,1)], "#CCCCCC")
  names(color_palette) <- status_levels
  
  # 构建可视化对象
  ggplot(rate_map_data) +
    geom_sf(aes(fill = status)) +
    # 颜色映射规范
    scale_fill_manual(
      values = color_palette,
      drop = FALSE,  # 强制显示所有类别
      na.translate = FALSE  # 已自行处理NA值
    ) +
    # 地图装饰元素
    annotation_north_arrow(
      location = "tl",
      style = north_arrow_fancy_orienteering(text_size = 8)
    ) +
    annotation_scale(
      location = "bl",
      width_hint = 0.3,
      style = "ticks"
    ) +
    # 图例配置
    guides(fill = guide_legend(
      title.position = "top",
      label.position = "bottom",
      nrow = 1,
      override.aes = list(color = "black", size = 0.5)
    )) +
    # 标签和主题
    labs(title = title) +
    theme_void() +
    theme(
      legend.position = "bottom",
      legend.title = element_blank(),
      plot.title = element_text(hjust = 0.5, face = "bold")
    )
}

# 打印所有地图（可选）
print(plots)

# 3.3.2 统计逐年状态数 
# 计算每年的状态数量
status_counts <- full_grid %>%
  group_by(year, status) %>%
  summarise(count = n(), .groups = 'drop') %>%
  complete(year, status, fill = list(count = 0))

status_counts <- status_counts %>%
  group_by(year) %>%
  arrange(desc(status)) %>% # 按照新设定的因子顺序排序
  mutate(cumulative_count = cumsum(count),
         label_position = lag(cumulative_count, default = 0) + count / 2) # 标签位于中间
writexl::write_xlsx(status_counts, "./output/表格5_status_counts.xlsx")

if (!"感染状态_计数" %in% names(wb2)) {addWorksheet(wb2, "感染状态_计数")}
writeData(wb2, sheet = "感染状态_计数", x = status_counts)

# 绘制堆叠柱状图并添加所有状态的标签
ggplot(status_counts, aes(x = factor(year), y = count, fill = status)) +
  geom_bar(stat = "identity", position = "stack") +
  geom_text(aes(y = label_position, label = count),
            vjust = 0.5, color = "black", size = 3, hjust = 0.5) + # 调整vjust使标签居中
  scale_fill_manual(values = c("从来没有感染过" = "lightgreen",
                               "曾经感染过" = "orange",
                               "当年首次感染" = "red",
                               "无病例报道" = "gray80"),
                    na.value = "gray70") +
  labs(title = "Infection Status of Counties Over Years",
       x = "Year",
       y = "Number of Counties",
       fill = "Status") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # 旋转x轴标签以提高可读性
##ggsave("./figs/图3.3.2 逐年状态区县数分布.png")

# 5.2 分阶段新发感染区县------
# 
# # 定义时间段
# time_periods <- list(
#   all_years = 2004:2023,
#   y2005_2015 = 2005:2015,
#   y2016_2020 = 2016:2020,
#   y2021_2023 = 2021:2023
# )

library(dplyr)
library(ggplot2)
library(sf)
library(writexl)
# library(ggsave)
library(stringr)

# 初始化函数以更新状态列
update_status <- function(df) {
  df %>%
    group_by(county_code) %>%
    arrange(year) %>%
    mutate(
      cumulative_infections = cumsum(!is.na(cases) & cases > 0),
      status = case_when(
        !is.na(cases) & cases > 0 & cumulative_infections == 1 ~ "当年首次感染",
        !is.na(cases) & cases > 0 ~ "曾经感染过",
        TRUE ~ "无病例报道" # 使用TRUE来捕获所有其他情况
      ),
      status = ifelse(is.na(cumulative_infections), "从来没有感染过", status)
    ) %>%
    ungroup() %>%
    select(-cumulative_infections)
}

# 获取所有年份和区县的组合，并初始化状态列
full_grid <- expand.grid(
  year = seq(min(yearly_incidence_county$year), max(yearly_incidence_county$year)),
  county_code = unique(yearly_incidence_county$county_code)
) %>%
  left_join(yearly_incidence_county, by = c("year", "county_code")) %>%
  mutate(status = "无病例报道") %>%
  update_status()

# 保存 full_grid 数据
write_xlsx(full_grid, "./output/表格5.1_感染状态_历年区县.xlsx")

# 绘制逐年地图函数
plot_status_map <- function(rate_data, title, map_data) {
  rate_map_data <- map_data %>%
    left_join(rate_data, by = c("GB1999" = "county_code"))
  
  centroids <- st_centroid(st_geometry(rate_map_data))
  rate_map_data$centroid_x <- st_coordinates(centroids)[,1]
  rate_map_data$centroid_y <- st_coordinates(centroids)[,2]
  
  ggplot() +
    geom_sf(data = rate_map_data, aes(fill = status)) +
    scale_fill_manual(values = c("从来不曾感染过" = "lightgreen",
                                 "曾经感染过" = "orange",
                                 "当年首次感染" = "red",
                                 "无病例报道" = "white"),
                      na.value = "gray50") +
    labs(title = title, fill = "Status") +
    theme_void() + 
    theme(
      plot.background = element_rect(fill = "white"),  
      plot.title = element_text(hjust = 0.5),          
      legend.position = "bottom",                      
      legend.box.background = element_rect(fill = "white", color = NA)  
    ) +
    annotation_north_arrow(location = "tl", which_north = "true",
                           style = north_arrow_fancy_orienteering) + 
    annotation_scale(location = "bl") + 
    geom_text(data = rate_map_data, aes(x = centroid_x, y = centroid_y, label = NAME), size = 3, vjust = -0.5)
}

# 分别绘制不同年份的地图并保存图片
dir.create("figs/逐年状态地图", showWarnings = FALSE, recursive = TRUE)
plots <- lapply(unique(full_grid$year), function(yr) {
  year_data <- filter(full_grid, year == yr)
  
  p <- plot_status_map(year_data, title = paste("Infection Status for", yr), map_data = map_data)
  
  ##ggsave(filename = file.path("figs/逐年状态地图", paste0("infection_status_", yr, ".png")),
         # plot = p, width = 10, height = 8, dpi = 300)
  
  p
})

# 计算每年的状态数量
status_counts <- full_grid %>%
  count(year, status, name = "count") %>%
  complete(year, status, fill = list(count = 0)) %>%
  group_by(year) %>%
  arrange(desc(status)) %>%
  mutate(cumulative_count = cumsum(count),
         label_position = lag(cumulative_count, default = 0) + count / 2)

# 保存每年状态数量
# status_counts %>% View()
write_xlsx(status_counts, "./output/表格5.1_感染状态_新发感染区县数量.xlsx")

# 定义所有可能的状态
status_counts$status %>% unique()
status_levels <- c("从来不曾感染过", "无病例报道", "曾经感染过", "当年首次感染")

# 绘制逐年堆叠柱状图
ggplot(status_counts, aes(x = factor(year), y = count, fill = status)) +
  geom_bar(stat = "identity", position = "stack") +
  geom_text(aes(y = label_position, label = count), vjust = 0.5, color = "black", size = 3, hjust = 0.5) +
  scale_fill_manual(values = c("从来不曾感染过" = "lightgreen",
                               "曾经感染过" = "orange",
                               "当年首次感染" = "red",
                               "无病例报道" = "gray90"),
                    na.value = "gray50",
                    limits = status_levels) +
  labs(title = "Infection Status of Counties Over Years",
       x = "Year",
       y = "Number of Counties",
       fill = "Status") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) 
ggsave("./figs/图3.3.2 逐年状态区县数分布.png")

# 分析并绘制各时间段的地图
dir.create("figs/分阶段状态地图", showWarnings = FALSE, recursive = TRUE)

plots_by_period <- lapply(names(time_periods), function(period_name) {
  years_in_period <- time_periods[[period_name]]
  
  period_data <- full_grid %>%
    filter(year %in% years_in_period) %>%
    update_status() %>%
    distinct(county_code, .keep_all = TRUE)
  
  p <- plot_status_map(period_data, title = paste("Infection Status for", period_name), map_data = map_data)
  
  ##ggsave(filename = file.path("figs/分阶段状态地图", paste0("infection_status_", period_name, ".png")),
         # plot = p, width = 10, height = 8, dpi = 300)
  p
})

# 每个阶段的状态数量
status_counts_combined <- bind_rows(lapply(names(time_periods), function(period_name) {
  years_in_period <- time_periods[[period_name]]
  
  period_data <- full_grid %>%
    filter(year %in% years_in_period) %>%
    update_status() %>%
    distinct(county_code, .keep_all = TRUE) %>%
    count(status, name = "count") %>%
    complete(status, fill = list(count = 0)) %>%
    mutate(period = period_name)
}))

# 保存结果
status_counts_combined
write_xlsx(status_counts_combined, "./output/表格5_感染状态_分阶段.xlsx")


if (!"感染状态_计数_分阶段" %in% names(wb2)) {addWorksheet(wb2, "感染状态_计数_分阶段")}
writeData(wb2, sheet = "感染状态_计数_分阶段", x = status_counts_combined)

saveWorkbook(wb1, file = "./output/0_发病率.xlsx", overwrite = TRUE)
saveWorkbook(wb2, file = "./output/0_莫兰指数.xlsx", overwrite = TRUE)

# # 绘制分阶段堆叠柱状图
# ggplot(status_counts_combined, aes(x = factor(period), y = count, fill = status)) +
#   geom_bar(stat = "identity", position = "stack") +
#   geom_text(aes(y = count / 2, label = count), vjust = 0.5, color = "black", size = 3, hjust = 0.5) +
#   scale_fill_manual(values = c("从来不曾感染过" = "lightgreen",
#                                "曾经感染过" = "orange",
#                                "当年首次感染" = "red",
#                                "无病例报道" = "gray80"),
#                     na.value = "gray50",
#                     limits = status_levels) +
#   labs(title = "Infection Status of Counties Over Time Periods",
#        x = "Time Period",
#        y = "Number of Counties",
#        fill = "Status") +
#   theme_minimal() +
#   theme(axis.text.x = element_text(angle = 45, hjust = 1)) 
# 
##ggsave("./figs/图3.3.2 分阶段状态区县数分布.png")

